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We investigate the possibility that nonlinear gravitational effects influence the preheating era after 
inflation. Our work is based on numerical solutions of the inhomogeneous Einstein field equations, 
and is free of perturbative approximations. The one restriction we impose is to limit the inhomo- 
geneity to a single spatial direction. We compare our results to perturbative calculations and to 
solutions of the nonlinear field equations in a rigid (unperturbed) spacetime, in order to isolate 
gravitational phenomena. We consider two types of initial conditions: where only one mode of the 
field perturbation has a non-zero initial amplitude, and where all the modes begin with a non-zero 
amplitude. Here we focus on preheating following inflation driven by a scalar field with a quartic 
potential. We confirm the broad picture of preheating obtained from the nonlinear field equations 
in a rigid background, but gravitational effects have a measurable impact on the dynamics for both 
sets of initial data. The rigid spacetime results predict that the ampfitude of a single initially excited 
mode drops rapidly after resonance ends, whereas in the relativistic case the amplitude is roughly 
constant. With all modes initially excited, the longest modes in the simulation grow much more 
rapidly in the relativistic calculation than with a rigid background. However, we see no evidence 
for the sort of gravitational collapse associated with the formation of primordial black holes. The 
numerical codes described here are easily extended to more complicated resonant models, which we 
will examine in the future. 
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I. INTRODUCTION 

For nearly a decade it has been known that coher- 
ent processes at the end of inflation can drive explosive 
particle production ^ via parametric resonance, lead- 
ing to an era of preheating in the post-inflationary uni- 
verse. Earlier assumptions about the transfer of energy 
from the inflaton to "regular" matter have been over- 
turned, and the phenomenology of the post-inflationary 
universe is thus much richer than previously suspected. 
However, the corollary of such interesting behavior is un- 
derstanding preheating has been a long and arduous task 
I |l5|,|l8|jl^jl|,|l|-|8). In particular, the problem is in- 
trinsically nonlinear and takes place far from thermal 
equilibrium. Moreover, employing simplifying assump- 
tions that are often useful in other circumstances can 
suppress crucial features of the preheating era, and when 
they are discarded the resulting predictions can change 
dramatically. 

The simplest model to exhibit parametric resonance is 
chaotic inflation driven by the quartic potential. 



Vi^) = J0^ 



(1) 



The (j) particles themselves are produced parametrically, 
so there is no need to couple the inflaton to other bosons, 
which is a prerequisite for resonance when V((j)) — 
771^0^/2. By shifting to conformal time and rescaling 
the field, all explicit dependence on the expanding back- 
ground can be removed from the field's equation of mo- 
tion. Consequently, the analysis of parametric resonance 



for X(f>'^ is effectively performed in Minkowski space, and 
is thoroughly discussed by Greene et al. ||lj . The insta- 
bility diagram contains a single, narrow resonance band 
and, to first order, co-moving modes which are initially in 
resonance appear to remain there for ever. Thus stochas- 
tic resonance, associated with a given wavelength passing 
through multiple resonance bands, does not occur. How- 
ever, when the back reaction of the created particles on 
the field evolution is incorporated the resonant growth 
terminates |0,|ll[. 

The inflaton field, (p, is described by a homogeneous 
background part, 4'o{t), and an inhomogeneous part, 
S(j){t, x'), which is usually expressed in terms of its Fourier 
modes, Scfik- Modes with wavenumbers, k, inside the 
resonance band grow exponentially, while other modes 
have oscillatory solutions. Because of the amplification of 
some of the S(j)k , the "matter" develops a substantial in- 
homogeneous component. The background spacetime is 
usually assumed to be rigid, and therefore unperturbed. 

A full treatment of parametric resonance must include 
the inhomogeneous metric induced by the inhomoge- 
neous matter. Incorporating metric perturbations into 
the analysis of parametric reson ance has bee n the subject 
of considerable recent attention 
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bation theory cannot give a full description of resonance, 
since once the 6(l)k grow large, mode-mode couplings be- 
come significant. Consequently, a better understanding is 
obtained from a direct (numerical) solution of the nonlin- 
ear field equations. These equations are nonlinear due to 
the higher order terms in the potential, which are a nec- 
essary condition for resonance. However, in these Simula- 



tions the background spacetime is an unperturbed FRW 
universe, so gravitational effects are implicitly ignored. 
To capture both the full nonlinear physics of preheating 
and any gravitational effects that influence the dynam- 
ics, we solve the evolution equations for the fields and the 
spacetime background. Our only a priori restriction is to 
confine the inhomogcneity to a single spatial dimension, 
which reduces the differential equations we must solve to 
a 1 + 1 dimensional system. 

Finally, it has been argued that second order terms in 
cosmological perturbation theory contribute to the over- 
all dynamics via back reaction effects [^,|0| . In the con- 
text of preheating the one calculation which explicitly 
includes second order effects showed that they did not 
significantly affect the dynamics ^^Q The broader is- 
sue of the back reaction of perturbations on cosmological 
spacetimes can be addressed with our code, and we in- 
tend to pursue this topic in the future. 

In general, we can anticipate two possible outcomes. 
If gravitational effects do not alter the dynamics the ini- 
tial phase of the evolution of the metric perturbations 
will reproduce the results of gauge invariant perturba- 
tion theory |l^Jl^ , pl| -p3[ , and at later times, when the 
perturbative limit is no longer valid, we will recover the 
results of the lattice simulations of the nonlinear field 
equations in a rigid background. On the other hand, if 
gravitational effects do play an important role, we will see 
phenomena that are not predicted by other approaches. 

Two possible gravitational effects have been identi- 
fied by previous work. On small scales, the increased 
inhomogeneity induced by resonance may lead to local 
gravitational collapse, and the formation of primordial 
black holes |22|. Since they decay by emitting Hawking 
radiation, any such holes would provide an additional 
channel for thermalizing the post infiationary universe 
[ p4| . Many inflationary models do lead to the forma- 
tion of primordial black holes, provided the perturbation 
spectrum they produce has more power at short scales 
than at COBE scales jSgJS^]. However, if preheating cre- 
ates a large amount of inhomogeneity at small scales, it 
could lead to primordial black hole formation in models 
where they would not otherwise be expected. On scales 
much larger than the post-inflationary Hubble length, it 
has been proposed that nonlinear effects enhance met- 
ric perturbations pq | , potentially undermining the usual 
inflationary predictions for structure formation and the 
anisotropy of the microwave background. Consequently, 
during our analysis we are particularly interested in look- 
ing for evidence for both gravitational collapse at short 
length scales, and the enhancement of metric perturba- 
tions at long length scales. Finally, any changes in the 
initial resonance structure due to the inclusion of gravi- 



tational perturbations would be revealed by our simula- 
tions. 

This paper extends the approach described by us in 
p4| (see also [^,Q) in two ways. Firstly, we generalize 
the allowed range of initial data so that Scj) can be an ar- 
bitrary function of position on the initial slice. Because 
of this, we can analyze the preheating era after chaotic 
inflation driven by a quartic potential, equation (|l|) , with 
a realistic set of initial conditions. The further extension 
to models with more than one scalar fleld is straightfor- 
ward and poses no new technical challenges. However, as 
the resonance provided by A0^ potential is conflned to a 
narrow band, it provides an excellent laboratory for in- 
vestigating both the perturbative and non-perturbative 
regimes, before we turn our attention to more complex 
two-field models, or single field models with more com- 
plex potentials. 

In the following section, we assemble familiar results 
that describe the evolution of a homogeneous universe 
after the end of inflation. We then introduce the evolu- 
tion equations for small perturbations, and the nonlin- 
ear equations of motion for the inhomogeneous flelds in 
an unperturbed spacetime. In Section 3, we derive the 
specific form Einstein field equations for the inhomoge- 
neous case, describe how we fix the initial conditions; we 
then outline the numerical algorithm we use to evolve the 
equations. Section 4 describes the results of our simula- 
tions. We look at two different sets of initial data. In 
the first, we start with a field perturbation that consists 
of single Fourier mode, chosen to lie inside the resonance 
band. This choice facilitates comparison with perturba- 
tion theory, by minimizing the effect of mode-mode cou- 
plings. The second, more realistic, choice of initial data is 
to begin with a field perturbation where all modes have 
a non-zero amplitude. Finally, in Section 5 we discuss 
our results and the questions they raise. For both the 
simulations, we find noticeable differences between the 
approximate calculations and the results obtained from 
solving the Einstein field equations. However, these do 
not appear to be large enough to overturn the overall 
picture obtained from lattice simulations of the nonlin- 
ear field equations in a static background. 



II. EVOLUTION EQUATIONS 

A. Homogeneous Evolution 

For infiation in a flat, homogeneous universe, we have 
the general result 



*However, a very recent paper by Bassett ei al. takes issue 
with these conclusions. 
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where, as usual, a{t) is the scale factor, and overdots 
stand for differentiation with respect to the time, t. The 
gravitational coupling is k^ = SttG = "Hii/m^y For the 
special case of a quartic potential, we make the transfor- 
mation 



where 
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a{t) 



dt. 



X = 00, 



after which the constraint, equation (0), becomes 
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and X obeys 



X" + Ax' 



-X-0. 



(5) 



(6) 



(7) 



Primes refers to differentiation with respect to the con- 
formal time, T]. It is easy to show that the X'^' I '^ ^^id 
Xa!' I a terms in equations (g) and (|7|) respectively quickly 
become negligible after the end of inflation ^^ . 
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(11) 



The time dependent quantities in equation (|lO|) are all 
evaluated at an (arbitrary) moment when the kinetic en- 
ergy of the inflaton vanishes. Strictly speaking, Greene 
et al. [Q do not analyze equation (Q), since they assume 
that only the field is perturbed, and that there are no 
metric perturbations. However, the resonance structure 
is not significantly changed by the presence of metric per- 
turbations. 

Because of the scaling properties of the solution, the 
values of k which fall inside the resonance band are not 
time dependent. This greatly simplifies the structure of 
resonance in this model, since only a small number of 
modes experience any resonant amplification. While we 
are interested the general question of nonlinear gravita- 
tional effects during preheating, this property of X(f>^ the- 
ory makes it an excellent testbed for exploring nonlinear 
gravitational effects during resonance, and the tools we 
are developing here will later be applied to more compli- 
cated models. 



B. Linear Perturbation Theory 

The onset of resonance can be studied perturbatively. 
In many cases, the perturbative analysis of resonance is 
based on an expansion about the homogeneous solution 
for the field, and assumes that the spacetime background 
is unperturbed. However, including metric perturbations 
(which couple to perturbations in the matter) does not 
make the problem a great deal more difficult, and this is 
the approach we take here. Perturbation theory is most 
conveniently formulated using Mukhanov's variable, Q 
[0,n8 31 33], which is related to (the fc-th Fourier mode 
of) the gauge invariant perturbation $fe [pl[ 



by 



Qk^' 



^' 



^<J>., 



(8) 



where (5(/)}5'' is the fc-th Fourier mode of the gauge in- 
variant perturbation of 0, and $ is the gauge invariant 
metric perturbation | |3l[ |. With n fields we have n dis- 
tinct Q, which obey the coupled, second order equations 
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= 0, (9) 



where the subscript on Q now labels the field, and not 
the Fourier mode. 

The only relevant resonance band is O] 



k^ = c^Xa^ 



1 0=0 



(10) 



C. Nonlinear Field Equations 

The next level of sophistication is to include all terms 
which are nonlinear in the fields, but to assume that the 
background spacetime is unperturbed. In the context of 
Xcj)* theory, this problem was first considered by Khleb- 
nikov and Tkachev M . The perturbative analysis of pre- 
heating involves ordinary differential equations, but the 
nonlinear field equations are partial differential equations 
which are typically solved numerically on a spatial grid. 
Resonance causes individual modes to become large, so 
mode-mode couplings are significant and first order per- 
turbation theory breaks down. The only nonlinear cou- 
plings ignored by these solutions are those attributable to 
gravity, and they thus provide a much more complete de- 
scription of resonance and preheating than perturbation 
theory. The purpose of this paper is to find out whether 
this description is complete, and we do this by compar- 
ing solutions of the nonlinear field equations with the full 
relativistic calculation we describe in the next section. 

We derive the nonlinear field equations by adding spa- 
tial gradient terms to equation (1^), 



X" - V^x + Ax' X = 0, 

a 



(12) 



where V = 9i and a prime denotes differentiation with 
respect to conformal time. Working with the conformal 
field, Xi simplifies the equations. This equation includes 
gradients of x, but the metric is implicitly assumed to be 
a FRW spacetime. However, if the spatial gradients are 
non-zero, the metric cannot be perfectly homogeneous. 



so gravitational effects on preheating are ignored when 
this equation is solved. 

Ignoring perturbations, the post-inflationary expan- 
sion mimics a radiation dominated universe, since on 
average the pressure and density of the oscillating field 
obey p ~ pi'i. In this limit, a(r\) ex 77, and the a" j a term 
rapidly becomes negligible after the oscillatory phase be- 
gins. 

Our numerical code solves both the Einstein field equa- 
tions and equation ([l2|). All the nonlinear terms in equa- 
tion (n2h also appear in the relativistic equations. By 
comparing x with the field evolution obtained from the 
relativistic calculation we can isolate gravitational effects 
during the preheating era. 

To recover from our solution for equation (|lj) we 
divide x by 0,. In order to obtain a, one could solve the 
ordinary differential equations which give a, or assume 
an exact radiation dominated expansion. In practice we 
extracted a by averaging the metric functions derived 
from our numerical solutions of the full Einstein equa- 
tions. Since the average of the perturbation vanishes, we 
recover the background solution in a self-consistent man- 
ner. However, we followed [Q and dropped the a" j a term 
from equation (|l^). 

III. EINSTEIN FIELD EQUATIONS 

We now turn our attention to the explicit form of the 
Einstein field equations which we will solve numerically. 
To reduce the computational complexity of the problem, 
we assume that the universe has a planar symmetry; the 
metric functions depend only on t and z, and are indepen- 
dent of X and y. In |2J], we examined possible nonlinear 
gravitational effects after rr?(^'^ inflation using the metric 

which describes an inhomogeneous universe in which the 
dx^ -I- dy^ sections have zero spatial curvature. In prin- 
ciple, we could retain this metric here, but in practice it 
is simpler to work with a metric for which the field equa- 
tions reduce to the conformal time version of the homo- 
geneous system, rather than physical time which is the 
homogeneous limit of equation (03) . The oscillation time 
for the homogeneous field is constant when expressed in 
conformal time, and this property persists in the nonlin- 
ear system. Thus the time-scale governing the evolution 
does not change significantly during the simulation, and 
we can use a constant timestep, simplifying our numerical 
problem. 

Specifically, a co-ordinate transformation from t and z 
to 77 and ^ allows us to re-write equation (nS) as: 

ds" = o?{j], C) (rf'y' - dC) - PHv, C) {dx^ + dy^). (14) 

With this metric, the non-trivial components of the Ein- 
stein tensor are 
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The field Lagrangian is, 
2 a2 



VW, 



and yields the following equation of motion for 
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The components of the energy-stress-momentum tensor 
are therefore 
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The corresponding field equations are 
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(18) 
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Finally, we have two constraints which must impose on 
the initial data, but which are then satisfied at all later 
times as a consequence of the field equations. 
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(21) 



A. Initial Conditions 

To solve the equations, we must specify the initial con- 
ditions. These cannot all be chosen arbitrarily, since the 
initial data must satisfy the constraint equations. We 
consider two possible scenarios. In the first we excite a 
single mode and track its evolution using the full non- 
linear equations. In the second, all modes are initially 
excited. Using the first approach it is easy to compare 
our results to a perturbative calculation. However, begin- 
ning with a configuration in which all the inhomogeneous 
modes have a small, non-zero amplitude is more realistic. 



1. A Single Excited Mode 

If we only wish to excite a single mode, we can pick ini- 
tial conditions which ensure that the constraints are au- 
tomatically satisfied, and thus do not need to be solved 
explicitly before beginning the actual integration. Fol- 
lowing the approach of M, we set a(0, Q — /3(0, Q — 1, 
c^(0, C) = 00- Equation (|21|) is solved if 



^r,(0,C) = 
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/3..(0,C) = l/y 



+ V{c^o) = C. 



(22) 



(23) 



Our choice of C ensures (a,,,(0, C)) = (/?,,, (0, C)), where 
(• • •) is a spatial average. The actual inhomogeneity is 
injected through the inflaton velocity, which we choose 
to be 



0,,(o,c) = n + 



/27rfcC\ 
e sin — - — . 



(24) 



where Z is the length of our "box" , and 11 is the average 
initial velocity. However, since we have assumed that 0, 
a and P are independent of C, on our initial slice we can 
only excite multiple modes if we assume that they are all 
initially in phase with each other. 



2. General Initial Data 

In the second scenario, where all the inhomogeneous 
modes are initially excited, we want to allow their ini- 
tial phases to be random. With just one mode the phase 
is irrelevant, but demanding that all modes begin with 
correlated phases would lead to an unacceptable loss of 
generality. In particular, we are interested in modeling 
the Gaussian fluctuations predicted by inflation, which 
are defined in part by having random phases. In this 
case there is no simple and obvious approach to solving 
the constraints. A general procedure does exist to fix 
the initial data in 3 + 1 dimensional numerical relativity 



P9[ , but it is not necessarily easy to implement. Since 
we are interested in a perturbation that is initially small, 
we will adopt the approach used by Goldwith and Pi- 
ran |4y], who simulate the onset of inflation in a 1 -I- 1 
dimensional universe. Their approach is to add a small 
amount of radiation, with a spatial dependence chosen 
to ensure that the overall density perturbation vanishes 
on the initial surface. Having done this the constraints 
become simple to solve, while the radiation introduced 
by this procedure always makes a trivial contribution to 
the overall density, and quickly decays away. 

We modify Goldwirth and Piran's approach slightly, 
and use an uncoupled scalar field, ip, instead of radiation. 
The energy density of a massless free field drops as a~^ in 
a homogeneous universe, in contrast to radiation which 
scales as a""*. This ensures that the contribution from ^ 
to the overall density drops more quickly than that of 0, 
which behaves like radiation. 

Including a second field in our code is simple. The field 
itself obeys the equation of motion. 
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(25) 



and we modify equations (|l8| ) and (|l9| ) by adding ^^^ 
and V".?) terms analogous to the derivatives of (p that are 
already there. Most importantly for us, the revised con- 
straints are 
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We solve these equations on the initial slice by choosing 
"0^^ and ip,ri so that the right hand sides of the constraints 
simplify, or 

l^l + l<Px + li^% + li^X + Va' = D\ (29) 

Now we are free to choose 
a(0, C) = PiO, = 1, a,^(0, C) = A^(0, C) = ^ (30) 



where H = Dk/^/3. Adding equation (27) to equa- 
tion (E9h gives two coupled equations for i/j^rj and ip^^, 
which we solve to find 
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(31) 
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where A^ = 20"^ - 2V - (0,,, + (j)^)'^ and B = -(j)x4'..r,, 
and we have used a(0,C) = 1- We can choose either of 
the two solutions. FinaUy, for these to be real, we need 
A^ — AB > 0, which is equivalent to demanding 



D > max 
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(33) 



We are free to choose D^, subject to the above constraint. 
In practice we want to minimize the impact of the extra 
field, so we set it close to the lower bound. However, our 
numerical scheme is somewhat more stable if the inequal- 
ity is not completely saturated, and in practice we chose 
a value of D somewhat larger than the minimum allow- 
able value. Finally, "0^^ and V',r; are periodic, as required 
by the boundary conditions. 



3. Choice of Initial Data 

We need to specify <j) and (/),,, on the initial slice. We 
begin our simulations at the end of inflation, so (j) is de- 
scribed by small fluctuations overlaid on top of the ho- 
mogeneous zero mode. The fluctuations are quantum in 
origin, and are generated by the usual inflationary mech- 
anism. In practice the spectrum will have some scale 
dependence, especially towards the end of inflation when 
the slow roll approximation breaks down. If we wished we 
could track the perturbations from their inception dur- 
ing inflation, but we will assume the spectrum is initially 
scale-free since our results do not depend on the precise 
form of the spectrum. 

Let q — (0(0, C)) and p = ((/)^^(0, C)), where (• • •) de- 
notes a spatial average; these values are chosen to be 
those at the end of inflation. We will write the perturba- 
tion in 0, for any time, as 



X = /3(0-(0)) 



(34) 



Since, the metric perturbations which are first order in 
77 vanish when 77 = 0, we can expand equation (1^) to 
obtain an equation for the homogeneous mode and a lin- 
ear equation for X, valid on the initial slice when ry = 0. 
Taking the Fourier transform of this latter equation gives 



^k,rm + (fc' +H^+ {Vc^c^ - K^V)) Xk = 0. 



(35) 



Thus for a short time near 77 = each fc-mode of X can 
be approximated as a simple harmonic oscillator with a 
fc-dependent frequency. 



l = k^ + H^ + {V^^-K^V). 



(36) 



Our approach will be to treat X{0,() as a free quan- 
tum field in the vacuum state; this amounts to the initial 
choice 
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(37) 
(38) 



where n is a number taken randomly from a Gaussian 
distribution of mean and variance 1, and r is taken 
randomly from the interval [0,1]. For equations ( |37| ) 
and (38) to be applicable the adiabatic limit must hold, 
namely ujk^r] '^ oj^- A more rigorous treatment is possi- 
ble but our approach is adequate for fixing the scale of 
the fluctuations. Moreover, we emphasize that while this 
choice of initial conditions is a good approximation when 
applied prior to the onset of resonance, its underlying 
assumptions will be strongly violated at later times. 

We insert </> into equation (36) to determine wj,, from 
which we learn X^ and, ultimately, 0. If we insist the 
fluctuations are smalffl we can solve this problem itera- 
tively. We approximate {V^^ — k'^V) by V^^{q) — K'^V{q). 
The final complication is that H also depends on 
through equation (Gfl) and equation (pOf). Our approach 
here is as follows: H^ — D^k^/3, and D^ = P4, + P4, is 
a constant, so D^ = (p^ -I- p^p). We approximate (p^) by 
p^ /2 + V{q) and then choose an e and insist (p^) — e(p0). 
Setting e ^ 1 ensures the initial energy density of ijj is 
small and fixes D and H . However, if e is too small it 
will not satisfy equation (p3|). This can only be checked 
when (j) and 0^^ have been determined: if e is too small 
and one must repeat the process with a larger value. In 
practice e can be on the order of ^, and the contribution 
of Tp rapidly becomes negligible as the universe expands. 

Once X{QX) and Xj^{Q, C,) have been found though, it 
is a simple matter to obtain the initial conditions for 
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and (f)^ri using equation (p4): 
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B. Coordinate Transformation 

To transform between the "conformal" and "physi- 
cal" co-ordinates, we need t = t(r?,C) and z = z{ri,()- 
The transformation that connects equation (n4) to equa- 
tion (13) reduces to the following pair of differential equa- 
tions 



t,ri = ^t'^^ + a'^ir],Q, 
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(41) 
(42) 



from which we find 



Hn a numerical scheme, requiring small fluctuations leads 
to a constraint on the lattice spacing. This is because the 
lattice provides a natural UV (and IR) cutoff to the unrenor- 
malized quantum field theory we have effectively introduced 
in equation ( pTI ) and equation (Pq). 
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and the specification is completed with the foUowing set 
of initial conditions, 



t(o,c) = o, z(o,c) = C 



(45) 



Since both t and z are transformed, slices of constant rj 
in the conformal frame are not mapped directly to slices 
of constant t in the physical frame. 

As we have already noted, the {rj, () co-ordinates are 
convenient because the typical timescale for the field os- 
cillations is roughly constant in this frame. Moreover, the 
degree of spatial variation is a{r], () is significantly less 
than in the corresponding A{t, z), which is enhances our 
numerical stability. In practice, we solve equations (Es) 
and (H) alongside with the field equations themselves, 
thereby obtaining solutions for both forms of the metric. 



C. Numerical Code 

The actual numerical code is essentially the same as 
that used in p4| . The simulations are evolved forwards in 
time on an N point spatial grid, with periodic boundary 
conditions. We use a fourth-order differencing scheme to 
express equations (18) and ( p^ ) as N ordinary differential 
equations, and solve these using a fourth-order Runge- 
Kutta integrator. Numerical noise is controlled with a de- 
aliasing algorithm, which requires filtering Fourier com- 
ponents with wavelengths less than 8 grid spacings at 
each timestep, so our actual spatial resolution is effec- 
tively 16 times larger than the separation of the individ- 



ual points, (5C ^l 42 1 . We track the numerical accuracy of 



the solutions by studying how closely the constraint equa- 
tions are satisfied, and have also compared solutions ob- 
tained directly from the metric equation (p^), and those 
obtained with equation ( |l4[ ) and the co-ordinate trans- 
formation equations (E^) and (pi). 

The Fourier transforms needed for de-aliasing are car- 
ried out using the Fast Fourier Transform algorithm. 
While this works for any value of A^, it is most efficient 
when N is an integer power of 2. 

We choose our units so that k^ = A = 1. Setting k^ 
defines the energy scale, but on first sight choosing A = 1 
may seem to result in a loss of generality. However, all 
terms in the equations of motion which do not involve V 
or V^ contain two derivatives with respect to rj and/or 
^ and the rescaling rj -^ rj/y^^ C ~^ C/V^ ensures that 
all the terms in the equations of motion are proportional 
to A. Consequently, our choice of A is equivalent to a 
rescaling of time and length. All the values of <^ and f] 
given in our results reflect this scaling. Finally, when 
we plot the spectra for variables such as $ and Q we 
normalize k so that fc = 1 corresponds to the mode whose 
wavelength is equal to the initial Hubble radius. 



D. The Perturbative Limit 

Consider the connection between equation (^^ and the 
familiar gauge invariant theory of cosmological pertur- 
bations [pl| . In conformal time, the most general metric 
that describes scalar perturbations to a FRW metric with 
flat spatial sections and whose spatial dependence is re- 
stricted to the ^ direction is 

ds^ - «(??)'{ (1 + 2^(r/, 0)drj^ - 2b{rj, O^dx'drj 

- [(1 - 2^^(77, 0)6,, + 2E{r,, C),^J]}dx'dx^ . (46) 

By comparing this with equation ( |l4|) and usingthe def- 
inition of the gauge invariant perturbation, $, |3^j2J] we 
obtain 



^,CC = "^.CC 



1 d 

a drj 



(EMCa), 



1 d 
a drj 



a^rjCt — P^nP 



H{a^ - (?) 



(47) 



(48) 



In order to obtain Q, we also need the gauge invariant 
perturbation, 54>'-'^^K For our metric pl|. 



^^''" 



(49) 



where 5(j) is the gauge dependent perturbation measured 
at constant rj. After differentiating twice with respect to 
(^, and again comparing equation (Bq) with equation (n4) 
we find that 



^^ 



'XC 



1 d0O r 

a^ drj 



/3A.-^(«' -/?')). 



(50) 



Our numerical simulations thus give the second deriva- 
tives of <I> and (5(/)''''' with respect to C,. To compare our re- 
sults with perturbation theory, we need the Fourier trans- 
forms of $ and (5(5!)'*'' . We find these by transforming <I>_^^ 
and 54>^fl, and multiplying by — fc^ which, thanks to the 
general properties of the Fourier transform, yields (E>j. and 
5(t)^^'\ We obtain the background variables (0(77), (l)o{v)) 
by averaging over the spatial grid at each value of rj. 
Once this is done, we have all the raw materials we need 
to calculate Qfc from equation (@), and make the compari- 
son with perturbation theory. In models with more than 
one field, each field has a gauge invariant perturbation 
analogous to (50<*'' . 



IV. RESULTS 

A. A Single Excited Mode 

We begin our analysis with an initial field field pertur- 
bation that consists of a single mode, chosen to lie within 




FIG. 1. The evolution of a single resonant mode is plot- 
ted as a function of rj. The top panel shows the evolution of 
Q^ predicted by perturbation theory. The second plot shows 
the evolution of Q^ derived from the full Einstein field equa- 
tions. The third panel superimposes the perturbative and 
non-perturbative results at the moment when perturbation 
theory breaks down. 
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FIG. 2. We plot the evolution of Q^ (x/*^) derived from 
the nonlinear field equations with a static background. The 
specific mode plotted is the same as that shown in fig. (mj, 
but the normalization is different since we dropped the a" /a 
term before solving equation (|l2|). 



the resonance band. The parameters for the integration 
are 

TV = 65336 == 216, St] = SC = 0.0722, hes = 583, 
(j)o = -0.58132, 00 = 

where kres denotes the mode number of the single excited 
mode, and fc = 1 corresponds to the mode whose wave- 
length is equal in size to the box. The values for (po and (f>o 
correspond to the first turn-around point after inflation. 
In the plots we normalize k so that fc = 1 corresponds 
to a mode with a wavelength equal to the initial Hubble 
horizon size, and fc^es ~ 1-2. The initial field values corre- 
spond to the field being at the end of the first oscillation 
it makes after the end of inflation. With this value of 
N, approximately 4000 distinct Fourier modes can be re- 
solved after de-aliasing. The spacing between the points, 
S(^, ensures that the simulation volume initially contains 
approximately 500 distinct Hubble volumes. As the co- 
moving wavelength of the resonant modes is slightly less 
than the Hubble length at the end of the inflation we can 
resolve modes with frequencies up to 7 times higher than 
the frequencies of modes in the resonance band. 

Fig. p) contrasts the perturbative evolution of our sin- 
gle excited mode, found by solving equation (^ , with the 
Qij, extracted from the relativistic calculation. Initially, 
perturbation theory is valid and Q increases exponen- 
tially. While perturbation theory predicts that exponen- 
tial growth continues indefinitely, non-perturbative ef- 
fects eventually terminate the resonance, as shown in the 
plots. The termination of resonance is not attributable 
to gravitational effects, and is well known from studies 
of equation ([l^).Fig. (Q) shows Q^ obtained by solving 
equation (|lj), which assumes that $ = and so Q re- 
duces to dcj). While the resonance terminates, the subse- 
quent evolution differs from the relativistic result, since 
in the absence of gravity the amplitude of the pertur- 
bation decreases after the end of resonance. Conversely, 
the Q^ obtained from the full relativistic calculation re- 
mains comparatively large after the resonant growth ter- 
minates. Presumably this is due to the self-gravity of 
the perturbation, which only enters perturbation theory 
at second order. 

We do not expect the solutions for x/a and ((> to overlap 
exactly since we have dropped the a" /a term from equa- 
tion (p^. This term quickly becomes irrelevant, but its 
absence during the initial stages of the evolution means 
that a perfect match is not possible, even with the same 
initial conditions. Moreover, this calculation includes the 
"compensating" field ^, even though it is not strictly nec- 
essary with a single mode and this also changes the initial 
normalization of (j), relative to x/*^- 

In addition to terminating the resonance, nonlinear ef- 
fects induce mode-mode couplings which transfer power 
to modes outside the resonance band. This can be seen 
in fig. (0), where we plot the power spectrum of the gauge 
invariant perturbation, <&. As we specify the initial per- 
turbation in a specific gauge, the corresponding gauge 
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FIG. 3. The power spectrum, |"l>fc| of the gauge-invariant 
metric perturbation is plotted against k. A single resonant 
mode is initially excited, but since we specify the initial per- 
turbation in a specific gauge, all modes of $ have some initial 
power. The plots show |$fc| at 77 = 0.22, 414.8 and 902.8. The 
a;-axis is scaled to that a mode with a wavelength equal to 
the initial Hubble radius has k = 1. 



invariant quantity actually possesses some power in all 
modes. As time progresses, the nonlinear couplings am- 
plify the "harmonics" of the resonant mode. This effect is 
also apparent when we solve equation (]l3) , but the quan- 
titative results differ from the general relativistic calcu- 
lation due to the observed decay of the perturbations in 
the absence of gravity. 



B. General Initial Conditions 

Having explored the nonlinear evolution found when a 
single mode of the perturbation was initially excited, we 
turn to the more realistic case in which the initial ampli- 
tude of all modes is non-zero. The initial amplitudes are 
set using equations ( |37| ) and (38), while the parameters 
for the simulation were 






FIG. 4. The power spectrum, |<I>fe|, of the gauge-invariant 
metric perturbation is plotted against k, with all modes ini- 
tially excited. The spectra are plotted at 77 =0.0183, 274.5, 
439.2, 750.3 and 915, respectively. 
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FIG. 5. The power spectrum of \Qk\, is plotted against k, 
with aU modes initiaUy excited. The spectra are plotted at 
the same times as in fig. m). 
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FIG. 6. The power spectrum of |xfc/i|i derived from a so- 
lution of equation (|l2|), is plotted against k, with all modes 
initially excited. The spectra are plotted at the same times 
as in fig. (^ and fig. (^ . 
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TV = 262144 = 2l^ 5t] ^ SC, ^ 0.0183. 
00 = -0.58132, 00 = 

The actual volume is the same as in the previous simula- 
tion, since the higher value of N is balanced by the lower 
value of 5C,, but the spatial resolution is improved. 

The amplitude of the initial perturbation here is sig- 
nificantly larger than the value obtained by extrapolat- 
ing the COBE normalization to the scales probed by our 
simulation. We made this choice to reduce the time re- 
quired for nonlinear effects to become important. How- 
ever, apart from prolonging the epoch during which the 
perturbative approximation is valid, repeating the cal- 
culation with an initial inhomogeneity that is 100 times 
smaller does not appear to make any significant difference 
to our results. 

The most important difference between the multi- 
mode case and the previous computation is in the be- 
havior of the longest wavelength modes. Recall that in 
the linear approximation (and in the absence of a tran- 
sient, decaying solution), and with k <C aH, 



cx 77, $ ex 1 / a{t')dt\ 



(51) 



so that Q is oscillatory. With a single excited mode 
nonlinear couplings amplified higher harmonics of the 
resonant mode, but the the long wavelength (low k) 
modes did not grow, and obey the relationships in equa- 
tion (pll). Conversely, with all the modes initially excited 
the longest modes undergo considerable growth, implying 
that the linearity assumption of equation ( |5l| ) is violated, 
and that mode-mode effects make a large contribution to 
the evolution of the longest modes. In fig. (Q) we plot 
|$fc| for the multi-mode case at five different times dur- 
ing the simulation. Initially, only resonant modes grow 
significantly. However, once nonlinear effects begin to 
transfer power to the harmonics of the modes in the res- 
onance band, we see that modes near k ~Q also start to 
grow. The resonant growth terminates, but the "bands" 
continue to broaden, and eventually merge to produce a 
continuous spectrum. Finally, at the end of resonance 
the Xfe/a modes do not decay, as they did when only a 
single mode was excited. 

We understand the evolution as follows. The initial 
growth of modes in the resonance band is explained by 
the perturbative analysis. The subsequent enhancement 
of the harmonics and the modes near fc = is due 
to mode-mode couplings. Second order terms transfer 
power from modes labeled by fci and k2 to the modes 
with fci -I- ^2 and |fci — fc2|. With just one excited mode, 
ki = k2, so we can produce the mode-doubling seeing 
in the previous simulation. However, for a single mode 
ki — k2 = 0, and modes with low k are not excited. Fi- 
nally, the broadening of the "bands" is (at minimum) a 
third-order effect. For instance, given three modes in one 
of the narrow bands, combinations like k = ki + k2 — k^, 



transfer power to modes adjacent to the band. This ac- 
counts for the observation that the higher harmonics are 
excited well before the bands themselves start to broaden. 

We also note that there is an alternative hypothesis 
which might explain the enhancement of the long wave- 
length modes. A conventional perturbative calculation 
predicts that these modes do not undergo resonance, but 
this calculation ignores the back reaction of the reso- 
nant perturbations on the oscillating homogeneous mode. 
This will induce a subtle change in the forcing term in the 
perturbation equations, which is provided by the oscil- 
lating field, and this correction could conceivably induce 
new resonances at long wavelengths. This is consistent 
with our results, since a single excited mode will have a 
much smaller back reaction than that produced by the 
full spectrum of excited modes, accounting for the ab- 
sence of this effect in the single mode simulation. More- 
over, the resonance structure of A^'* inflation is known to 
be very sensitive to small changes in the driving function. 
If this is correct, the harmonics of the resonant modes 
would still be excited through mode-mode couplings, but 
the longest modes would be excited by a co-operative ef- 
fect which would involve all the resonant modes. We 
stress that we have not yet attempted to verify this hy- 
pothesis. Note too that the possibility of backreaction 
by perturbations on the homogeneous mode is discussed 
(albeit in different models) in p7p8]. 

The most important question is the extent to which 
these effects are attributable to nonlinear gravitational 
couplings, as opposed to the nonlinear couplings pro- 
vided by the Xcj)'^ potential. Looking at the solutions 
of x/o found via equation (|lj), we see a similar, but not 
identical, resonance structure to that predicted by rela- 
tivistic calculation. The discrepancy is most obvious at 
long wavelengths. In fig. (||) we plot Qk and in fig. (|) 
we show Xfe/a at the same times as the plots in fig. (B). 

It is more appropriate to compare x/a to Q (or (50''''') 
than to (f>, since the latter vanishes if the metric is unper- 
turbed. At short wavelengths, the Qk and Xk/o. are very 
similar, but at the longest wavelengths Qk and (50**'' un- 
dergo much more growth than Xk/a. Moreover, while the 
longest modes of Xk/o. do grow, albeit less dramatically 
than the corresponding Qk, the growth begins signifi- 
cantly later. This is illustrated by fig. (R), which plots 
the evolution of the longest modes as a function of j]. 
Conversely, fig. (||) shows the same plot for A: = 50, and 
in this case the difference between the evolution of x/a 
and the relativistic quantities is much less pronounced. 
The final plot in fig. (H) shows the growth of Q for a 
mode inside the resonance band (this is actually the mode 
that was excited in the single mode case). Looking at 
this plot, we see that the near exponential growth of the 
long-wavelength modes ceases when the resonance is ter- 
minated. 

We complete our analysis of figs (R) and (H) by numer- 
ically estimating the slopes of the the exponential seg- 
ments of the plots. These slopes contain an arbitrary 
multiplicative factor, fixed by the normalization of a. 
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which enters via the definition of the conformal time, 
rj. Numerically, and with our normalization of rj, we 
find that the exponential segments of the first four plots 
in fig. (0) have a slope of approximately 0.036. Con- 
versely, the resonant mode shown in the bottom panel 
has a slope that is half of this value (0.017), supporting 
our contention that the growth in the longest modes is 
due to nonlinear couplings between two adjacent resonant 
modes. Finally, the exponential segment of the plots in 
fig. (^) has a slope of around 0.08, confirming that this 
growth is driven by a strongly nonlinear multi-mode cou- 
pling. 

Since we are working with the full Einstein field equa- 
tions we can examine non-perturbative measures of cur- 
vature. For instance, C^ = C^^^yxC^'^^^, where C is the 
usual Weyl tensor, is a covariant scalar density that is 
sensitive to the local spacetime curvature. Consequently, 
we can define 






(52) 
(53) 



where R is the usual curvature scalar. These integrals are 
not covariant, since they depend on a specific co-ordinate 
choice. However, C'^/TZ^ is dimensionless, and gives a 
qualitative measure of the inhomogeneity on our hyper- 
surfaces, which we plot in fig. (^, because C^ vanishes 
everywhere in an FRW spacetime. Thanks to our choice 
of initial conditions, the metric functions are constants at 
?7 = 0, so C^ vanishes on the initial hypersurface. Hence, 
the magnitude of C^/R^ is a measure of how far our 
spacetime departs from a perfectly FRW universe. The 
initial entropy perturbation quickly generates a metric 
perturbation. This remains constant until the nonlinear 
growth of the longest modes sets in at around r] = 200, 
whereupon C^ /Ti? grows exponentially until the under- 
lying resonance is halted. However, after resonance ends 
C^ I'R? remains roughly constant. 

If the enhanced density perturbations generated dur- 
ing resonance led to localized gravitational collapse and 
the formation of primordial black holes (or "black walls" 
in the case we consider here), this would be obvious from 
C^ /Ti? as the curvature would diverge in a collapsing re- 
gion. Since this does not happen in our simulation we 
find no evidence that preheating after Xcj)'^ inflation leads 
to primordial black hole formation. However, this result 
should be interpreted with care. After inflation, the equa- 
tion of state resembles a relativistic gas, where p = p/3, 
and as a consequence gravitational collapse is much less 
efficient than in dust filled universe with p = Ojj Con- 
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■'■Recall the absence of structure formation in a hot dark mat- 
ter dominated universe, as compared to a cold dark matter 
universe with the same primordial perturbation spectrum. 



FIG. 7. From top to bottom, these plots show l^il, |Qi|, 
'/'ill I'^Xi/'^l 3-rid IQsssl as functions of 77. 
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FIG. 8. From top to bottom, these plots show |"l>5ol, IQsol, 
|5(^5o|, and |(5x5o/t[| as function of r?. 

sequently, even if preheating does sometimes lead to pri- 
mordial black holes, it would not be surprising that they 
were absent from the A^"* case. In future work we will 
study more general potentials, and will be able to con- 
sider this topic in greater detail. 

Apart from the length scale provided by the resonance 
band, the other important scale in the post-inflationary 
universe is the Hubble scale, H~^. However, this length 
is not constant during the course of the simulation. In 
fig. (H) we plot the scale factor corresponding to the un- 
perturbed FRW limit of our simulation, and we see that 



FIG. 9. This plot shows the scale factor a{ri), obtained by 
averaging the metric functions obtained from our simulation 
to remove the perturbation. 
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FIG. 10. The growth of C^/7^^ during the course of the 
simulation is plotted as a function of rj. 



the universe grows 90 times larger during the simulation. 
During the same interval the Hubble radius grows lin- 
early from an initial value of 10 to a final value of 900 (in 
our units). Our simulation region has a comoving width 
of 4797 (2^^(50, so the fc-th mode has a comoving wave- 
length 4797/fc. Modes with k ^ 500 are initially outside 
the horizon. At the time when the fc = 1 (77 ~ 200) mode 
begins to grow, only modes with k ^ 25 are outside the 
Hubble horizon. By the end of the simulation, only four 
or five modes have wavelengths greater than H^^. 

Modes with fc > 50 show few differences between the 
evolution of Qj, obtained from the relativistic calcula- 
tion, and the analogous result for Xk/a. Consequently, on 
these scales (which can safely be termed "sub-Hubble" ) 
we see no evidence that nonlinear gravitational effects 
infiuence the overall progress of preheating after Xcf)"^ in- 
flation. 

For the longest modes there are clear differences be- 
tween the results found in a fixed FRW background and 
the full relativistic calculation. Both the relativistic cal- 
culation and the nonlinear field equations in an unper- 
turbed background predict that the longest modes of the 
field perturbation grow significantly. However, the rela- 
tivistic results predict that this growth begins earlier and 
is an order of magnitude larger than that found from the 
field evolution in a rigid background. Note that in both 
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FIG. 11. The evolution of 7, equation (p4), which quanti- 
fies the accuracy with which we satisfy the constraint, equa- 
tion (CT), is plotted for our numerical solution. For an exact 
solution 7 would be identically zero. 



cases we are seeing the growth of modes that he outside 
the instantaneous Hubble horizon. 

At present, we do not identify the source of this dis- 
crepancy, but several possibihties exist. The most inter- 
esting explanation would obviously be back reaction, or 
some other nonlinear effect ||2^,g^,^,^. However, we 
stress that we have no evidence that the specific results 
of our simulations are predicted by these authors. 

There are a number of less glamorous possible explana- 
tions for the enhanced inhomogeneity, which we briefly 
review. Firstly, it could be a side effect of using peri- 
odic boundary conditions, which induce an artificial long- 
range correlation into our simulation. Similarly, we as- 
sume that the inhomogeneity is limited to one spatial 
dimension, and the nonlinear effects we observe may be 
an artifact of this restriction. In many nonlinear sys- 
tems the dimensionality has a crucial influence on the 
behavior, so we would still need to verify that any cos- 
mologically significant changes to long wavelength modes 
persisted in the fully inhomogeneous case. On the other 
hand, our calculation of x/o is also restricted to a sin- 
gle inhomogeneous dimension, so the dimensional depen- 
dence of Q and x/a would have to be different if this 
explanation is correct. Likewise, the three dimensional 
simulations of Khlebnikov and Tkachev |0] , show a band 
structure that is qualitatively similar to that obtained 
here for x/^i which suggest that at shorter wavelengths 
our results are reliable. However, since we consider with 
a 1 -|- 1 dimensional system we are able to resolve a large 
range of length-scales, whereas the simulation region in 
Khlebnikov and Tkachev's work is not significantly larger 
than the initial Hubble volume. Consequently, we cannot 
compare our calculation to theirs at super- Hubble scales. 

We also need to be sure that the results we have seen 
do not depend on our introduction of the -0 field to sat- 
isfy the initial constraint. Ideally in future work we will 
solve the constraints directly on the initial surface and 
dispense with ^ entirely. However, we have varied the 
initial fraction of ij: and checked that changing it made 
no significant difference to our results. 



A different possibility is that our "reference" hyper- 
surfaces, defined as the 3-planes with constant 77, are 
degenerate in some way. If this is true the degeneracy 
must be shared by ^(77, Q and ^(77, C,), which we compute 
from equations (O) and (E2|), as z and t remain smooth 
functions of Q and 77 during the simulation. Moreover, 
we have been careful to work with variables which are 
invariant (to first order) under gauge transformations. 

Numerically, we checked our solutions by ensuring that 
they did not change significantly when we altered the pa- 
rameters that govern the numerical integration scheme. 
In particular, the results appear to be independent of 
N , the number of points on our grid, the distance be- 
tween them, 5Q, and the timestep 5rj. The accuracy with 
which we fit the constraints provides an independent test 
of our results, and the constraints are indeed well satis- 
fied, especially during the era when nonlinear effects first 
become apparent. 

The constraints, equations (|2^) and (p7|), allow us to 
form an independent measure of the accuracy of our code. 
To measure how well our solution obeys a given con- 
straint we evaluate 



1' = 



JdCil 



./JdCP./Jdc'r 



(54) 



We have two different 7, corresponding to the two con- 
straints, while I and r are the left and right hand sides of 
the equation equation (Effl or equation ( p7| ) (depending 
on which constraint we are testing), and 7 is the aver- 
age fractional difference between the left and right sides 
of the constraint equations. In fig. (^l]) we plot the 7 
derived from equation (27), and we see that it is always 
small. Moreover, 7 decreases with Srj and 5(1. In general, 
the value of 7 derived from equation (|2^) is smaller than 
that obtained from equation (pTf), and we do not show it 
here. 

Our choice of N and SC must satisfy two opposing cri- 
teria. Ideally, we would like to be able to resolve modes 
which were much larger than the Hubble radius for the 
entire duration of resonance, but to prevent the build up 
of numerical error we must resolve all the harmonics of 
the resonant modes that are excited by the nonlinear cou- 
plings. This sets an upper limit on d(^ and Srj. Including 
super-Hubble modes sets a lower limit on the total size of 
the box. The combination of these two effects is to ensure 
that N must be quite large. Consequently, while it is fea- 
sible to include modes which are a few times larger than 
the Hubble radius at the end of resonance, looking at 
modes 100 times larger (or more) would require a signif- 
icant amount of computer time. A better approach may 
be to look for similar effects in models where the reso- 
nance is stronger, which means that the overall growth of 
the universe during the resonant era will be much smaller 
than in the A^** case. 
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V. CONCLUSIONS 
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We have analyzed parametric resonance after X(f)* infla- 
tion for two distinct sets of initial conditions. In the first 
we gave only one Fourier mode of the perturbation a non- 
zero initial amplitude, and in the second all modes were 
initially excited. The overall picture of the resonant era 
derived from our calculations, which incorporates non- 
linear gravitational effects, is similar to that obtained 
from analyses of the nonlinear field dynamics which as- 
sume that the background spacetime is rigid. However, 
in both simulations we saw new phenomena which appear 
to be due to nonlinear gravitational effects. With a sin- 
gle mode excited, the perturbation calculated in a rigid 
spacetime background decayed after the resonant growth 
terminated, but its amplitude after resonance remained 
large when nonlinear gravitational effects were included. 
When all modes were initially excited power was trans- 
ferred to long wavelength modes by second order effects 
in both the rigid background and nonlinear gravitational 
calculations, but in the latter case the overall growth is 
larger and begins sooner. However, we do not see any 
evidence for the formation of primordial black holes and 
the long wavelength modes cease growing once the expo- 
nential increase of the modes inside the resonance band 
comes to a halt. 

Given the recent interest in nonlinear gravitational ef- 
fects during preheating and in gravitational back reac- 
tion, the enhanced inhomogeneity we see at long wave- 
lengths when all modes are initially excited is a tantaliz- 
ing result. However, more work will be needed to confirm 
that these effects have a physical origin, and are not an 
artifact of our numerical scheme or underlying assump- 
tions. Moreover, if the effect is real, its precise origin 
and physical consequences remain to be determined. In 
addition, the symmetry condition we have imposed in 
our simulation can be shown to preclude the existence of 
tensor perturbations to the metric, so we cannot reach 
any conclusions about the resonant production of gravi- 
tational waves after inflation. 

Finally, the numerical codes and analytic techniques 
developed during this investigation are readily applica- 
ble to other resonant models, including those based on 
two (or more) interacting scalar fields. Since the struc- 
ture of the resonance after A0^ inflation is very simple, 
we can form a clear picture of the evolution of the in- 
homogeneous parts of the field and the spacetime metric 
during the resonant era. Consequently, in addition to its 
intrinsic interest, A^** is an excellent model to begin with 
as it has allowed us to perfect the techniques we will use 
in the analysis of other, more complicated, systems. This 
work is currently in progress, and its results will shed fur- 
ther light on the potential nonlinear gravitational effects 
we have uncovered here. 
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